function                                solution = solve_idsds(pars,Nu,loud)

% Reduced-form parameters        
pars.A_M                                = 1/(pars.r+pars.k+pars.theta);
pars.A_X                                = 1/(pars.r+2*pars.k)*pars.C;
pars.A_N                                = pars.k*pars.A_X;
pars.uPhi0                              = pars.M_c -(pars.r+pars.k+pars.theta)/(pars.r+pars.k)*(pars.M_c - pars.M_e);           
pars.oPhi1                              = pars.uPhi0 + pars.A_X/pars.A_M*( 1  + (pars.k)/(pars.r+pars.k) );    

if loud
fprintf('Model parameters \n');
fprintf('r                              = %5.3f \n',pars.r);
fprintf('k                              = %5.3f \n',pars.k);
fprintf('M_e                            = %5.3f \n',pars.M_e);
fprintf('M_c                            = %5.3f \n',pars.M_c);
fprintf('sigma                          = %5.3f \n',pars.sigma);
fprintf('theta                          = %5.3f \n',pars.theta);
fprintf('\n');       
end;
  
% Discretization parameters
Nu.X                                    = linspace(0,1,Nu.NX);                                                                                          % Grid for X
Nu.hX                                   = min(diff(Nu.X));                                                                                              % Step size on the grid for X
Nu.NM                                   = 2*Nu.N+1;                                                                                                     % Number of gridpoints for M
Nu.M_min                                = pars.M_c - sqrt(Nu.N/Nu.theta_target)*pars.sigma;                                                             % Lower bound for M grid
Nu.M_max                                = pars.M_c + sqrt(Nu.N/Nu.theta_target)*pars.sigma;                                                             % Upper boudn for M grid
Nu.M                                    = linspace(Nu.M_min,Nu.M_max,Nu.NM)';                                                                           % Grid for M
Nu.hM                                   = pars.sigma/(sqrt(Nu.N*Nu.theta_target));                                                                      % Step size on the grid for M
Nu.Xgrid                                = repmat(Nu.X,[Nu.NM,1]);                                                                                       % Meshed grid for X
Nu.Mgrid                                = repmat(Nu.M,[1,Nu.NX]);                                                                                       % Meshed grid for M
Nu.Deltat                               = 1/(Nu.N*Nu.theta_target);                                                                                     % Time step, fraction of a month
Nu.tol                                  = (pars.r+pars.k)*Nu.Deltat*1e-4;

if loud
fprintf('Discretization parameters \n');
fprintf('theta_target                   = %5.3f \n',Nu.theta_target);
fprintf('NX                             = %5.0f \n',Nu.NX);
fprintf('hX                             = %5.3f \n',Nu.hX);
fprintf('N                              = %5.0f \n',Nu.N);
fprintf('N_M                            = %5.0f \n',Nu.NM);
fprintf('M_min                          = %5.3f \n',Nu.M(1));
fprintf('M_max                          = %5.3f \n',Nu.M(end));
fprintf('uPhi(0)                        = %5.3f \n',pars.uPhi0);
fprintf('oPhi(1)                        = %5.3f \n',pars.oPhi1);
fprintf('dt                             = %5.3f months, = %5.3f days \n',Nu.Deltat,Nu.Deltat*30);
fprintf('\n');
end;

if loud
fprintf('Condition on X grid size \n');
fprintf('min(1/(r+2k),1/(N_X*k))        = %5.3f \n',min(1/(pars.r+2*pars.k),1/(Nu.NX*pars.k)));
fprintf('dt                             = %5.3f \n',Nu.Deltat);
if Nu.Deltat < min(1/(pars.r+2*pars.k),1/(Nu.NX*pars.k))
        fprintf('Condition on X grid size satisfied \n');
else
        fprintf('Condition on X grid size not satisfied \n');
        error;
end;
fprintf('\n');
else
if Nu.Deltat > min(1/(pars.r+2*pars.k),1/(Nu.NX*pars.k))
        Nu.Deltat
        1/(pars.r+2*pars.k)
        1/(Nu.NX*pars.k)
        error;
end;
end;

if loud
fprintf('Conditions on M grid size \n');
fprintf('1/(N_M*theta_target)           = %5.3f \n',1/(Nu.NM*Nu.theta_target));
fprintf('dt                             = %5.3f \n',Nu.Deltat);
if Nu.Deltat <= 1/(Nu.N*Nu.theta_target) & Nu.M_min < pars.uPhi0 & Nu.M_max > pars.oPhi1
        fprintf('Conditions on M grid satisfied \n');
else
        fprintf('Conditions on M grid not satisfied \n');
        error;
end;
fprintf('\n');
else
if ~( Nu.Deltat <= 1/(Nu.N*Nu.theta_target) & Nu.M_min < pars.uPhi0 & Nu.M_max > pars.oPhi1 )
        error;
end;
end;

%%%%%%%%%%%%%%%%%
% Solution setup
%%%%%%%%%%%%%%%%%

% Transition matrices for X
dX_m                                    = -pars.k*Nu.Deltat*Nu.X;
J_m                                     = diag(ones(Nu.NX,1),0) + (Nu.hX)^(-1)*diag(dX_m,0) + (Nu.hX)^(-1)*diag(-dX_m(2:end),1);
dX_p                                    = pars.k*Nu.Deltat*(1-Nu.X);
J_p                                     = diag(ones(Nu.NX,1),0) + (Nu.hX)^(-1)*diag(-dX_p,0) + (Nu.hX)^(-1)*diag(dX_p(1:(end-1)),-1);

% Transition matrices for M
u                                       = (pars.sigma^2 + Nu.hM*pars.theta*(pars.M_c - Nu.M))/(2*Nu.hM^2)*Nu.Deltat;
d                                       = (pars.sigma^2 - Nu.hM*pars.theta*(pars.M_c - Nu.M))/(2*Nu.hM^2)*Nu.Deltat; 
if pars.theta == 0
   u(1)                                        = 0;
   d(1)                                        = 0;
   u(end)                                      = 0;
   d(end)                                      = 0;     
end
J_M                                     = diag(u(1:(end-1)),1) + diag(d(2:end),-1) + diag(1-u-d,0);

% Flow benefits matrix
Flow                                    = (pars.M_e - Nu.Mgrid + pars.C*Nu.Xgrid)*Nu.Deltat;

tt = tic;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Iterations from below
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Solution under no adoption
B_worse                                 = 1/(pars.r+pars.k+pars.theta)*(pars.M_c - Nu.Mgrid) - 1/(pars.r+pars.k)*(pars.M_c - pars.M_e) ...
                                                + pars.C/(pars.r+2*pars.k)*( Nu.Xgrid + 0*pars.k/(pars.r+pars.k) );
Phi_worse                               = pars.uPhi0 + pars.A_X/pars.A_M*Nu.X;

% Iterate, starting from below
B_initial                               = B_worse;
crit                                    = 1;
n                                       = 0;
B_prev                                  = B_initial;
while crit > Nu.tol   
        Ap                                      = B_prev >= 0;
        Am                                      = ones(Nu.NM,Nu.NX) - Ap;
        B                                       =  Flow + (1-(pars.r+pars.k)*Nu.Deltat)*( Ap.*(J_M*B_prev*J_p) + Am.*(J_M*B_prev*J_m) );
        n                                       = n+1;
        crit                                    = max(max(abs(B-B_prev)));
        B_prev                                  = B;
end;
B_below                                 = B;
Phi_below                               = NaN(Nu.NX,1);
for nx=1:Nu.NX
        if B_below(1,nx) <= 0
                Phi_below(nx)                           = Nu.M(1);
                error('Min value of M is too high');
        else if B_below(end,nx) >= 0
                Phi_below(nx)                           = Nu.M(end);
                error('Max value of M is too low');
                else 
                        nm                                      = find( B_below(:,nx) < 0,1,'first' );
                        if isempty(nm) || (nm == 1)
                                error('Unclear error with the threshold');
                        else
                                Phi_below(nx)                           = Nu.M(nm-1) - Nu.hM/(B_below(nm,nx)-B_below(nm-1,nx))*B_below(nm-1,nx);
                        end;
                end;
        end;
end
solution.B_worse                        = B_worse;
solution.Phi_worse                      = Phi_worse;
solution.B_below                        = B_below;
solution.Phi_below                      = Phi_below;

%%%%%%%%%%%%%%%%%%%%%%%%
% Iterations from above
%%%%%%%%%%%%%%%%%%%%%%%%

% Solution under no adoption
B_best                                  = 1/(pars.r+pars.k+pars.theta)*(pars.M_c - Nu.Mgrid) - 1/(pars.r+pars.k)*(pars.M_c - pars.M_e) ...
                                                + pars.C/(pars.r+2*pars.k)*( Nu.Xgrid + 1*pars.k/(pars.r+pars.k) );
Phi_best                                = pars.uPhi0 + pars.A_X/pars.A_M*Nu.X + pars.A_X/pars.A_M*pars.k/(pars.r+pars.k);

% Iterate, starting from below
B_initial                               = B_best;
crit                                    = 1;
n                                       = 0;
B_prev                                  = B_initial;
while crit > Nu.tol   
        Ap                                      = B_prev >= 0;
        Am                                      = ones(Nu.NM,Nu.NX) - Ap;
        B                                       =  Flow + (1-(pars.r+pars.k)*Nu.Deltat)*( Ap.*(J_M*B_prev*J_p) + Am.*(J_M*B_prev*J_m) );
        n                                       = n+1;
        crit                                    = max(max(abs(B-B_prev)));
        B_prev                                  = B;
end;
B_above                                 = B;
Phi_above                               = NaN(Nu.NX,1);
for nx=1:Nu.NX
        if B_above(1,nx) <= 0
                Phi_above(nx)                           = Nu.M(1);
                error('Min value of M is too high');
        else if B_above(end,nx) >= 0
                Phi_above(nx)                           = Nu.M(end);
                error('Max value of M is too low');
                else 
                        nm                              = find( B_above(:,nx) < 0,1,'first' );
                        if isempty(nm) || (nm == 1)
                                error('Unclear error with the threshold');
                        else
                                Phi_above(nx)           = Nu.M(nm-1) - Nu.hM/(B_above(nm,nx)-B_above(nm-1,nx))*B_above(nm-1,nx);
                        end;
                end;
        end;
end
solution.B_best                         = B_best;
solution.Phi_best                       = Phi_best;
solution.B_above                        = B_above;
solution.Phi_above                      = Phi_above;

solution.pars                           = pars;
solution.Nu                             = Nu;

tt = toc(tt);
if loud
fprintf('ISDS time: %d minutes and %4.2f seconds\n', floor(tt/60), rem(tt,60));
end;

end
